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Abstract Time-like orbits in Schwarzschild space-time are presented and classi- 
fied in a very transparent and straightforward way into four types. The analytical 
solutions to orbit, time, and proper time equations are given for all orbit types in 
the form r = r(A), t = t(x), and t = t(x), where A is the true anomaly and \ IS a 
parameter along the orbit. A very simple relation between A and x is a l so shown. 
These solutions are very useful for modelling temporal evolution of transient phe- 
nomena near black holes since they are expressed with Jacobi elliptic functions 
and elliptic integrals, which can be calculated very efficiently and accurately. 
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1 Introduction 

When modelling physical phenomena occurring in strong gravitational field of 
black holes, it is common to work in the Schwarzschild space-time 1,2,3,4,5,6,7, 
8 9, 10 . The same applies also for determination of orbital parameters of objects 
orbiting so close to a black hole that the orbits are affected by relativistic effects, 
e.g. highly eccentric S stars at the Galactic Centre [111112] . For this purpose, an 
efficient and accurate way for calculating time-like orbits in the Schwarzschild 
space-time is required. Moreover, if we are interested in time-dependence of these 
phenomena, a method for solving the time equation is also required to calculate 
the temporal evolution and the dynamics. 

Since in such models, the number of calculations can rapidly increase either 
because of increasing the number of points, or extending the time of the simu- 
lation, or reducing the time-step size, it is very desirable to have a very efficient 
and accurate method for solving these equations. For example, in a model which 
includes time-dependant gravitational lensing, it is easy to miss the moment of 
the strongest lensing when the Einstein ring appears, if the time-step is too large. 
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Consequently, the calculated signal, as received by a distant observer, lacks this 
distinctive characteristic. 

Although the numerical integration or post-Newtonian approximation yield 
useful results in specific cases, analytical solutions of the orbit and time equation 
are simpler, more efficient regardless of the accuracy required (as shown by Delva 
[13]). and can be used in all cases (weak field limit, strong field limit). The well 
known work of Chandrasekhar [T3] and Rauch [15] , where the solutions to geodesic 
equations are expressed in terms of elliptic integrals, has been followed by Cadez 
16,17 and Gomboc [18] who inverted the expressions of Chandrasekhar [2] and 
Rauch [15] into Jacobi elliptic functions, which no longer contain the branch ambi- 
guity. For light-like geodesies, Cadez and Kostic [T7] presented a very simple and 
straightforward way of characterizing the orbits which depend only on one param- 
eter, as well as giving analytical solutions to the time equation and a method of 
determining a light-like geodesic between two arbitrary points (and thus facilitat- 
ing ray-tracing used in numerical modelling of dynamical phenomena near black 
holes) . 

Cruz et al. [19] have classified the light-like and time- like geodesies according to 
the effective potential and found similar analytical solutions as [17 , however they 
give solutions to time equation only for radial and circular orbits. Hioe and Kuebel 
|20] present analytical solutions to orbit equations, classify them according to two 
parameters, and show extensive tables of different values of these parameters for 
corresponding orbits. They, however, do not give any solution to time equation. 

To complement previous work on light-like orbits [T7] , the complete analytical 
solutions of the time-like geodesies and the time equation for all orbit types are 
presented in this paper: in the form r = r(A) (where A is the true anomaly) for the 
radial coordinate r, and in the form t = t(x) and r = t(x) f° r time t and proper 
time r, respectively, with a very simple relation between A and x- 



2 Schwarzschild space-time 

In Schwarzshild space-time we use Schwarzshild coordinates t, r, 6, ip. The Hamil- 
tonian, from which geodesic equations are derived is 



where p^ are canonical momenta and natural units c = G = 1 are used. The 
constants of motion are: value of Hamiltonian (H) and Lagrangian (L), energy 
E = pt, three components of angular momentum (1), longitude of periapsis (w), 
and time of periapsis passage (t p ). For time-like geodesies, the value of Hamiltonian 
is = —1/2. 

In order to describe the position along the orbit, as well as the orientation of 
the orbit, we introduce another local inertial (right-handed) orthonormal tetrad 
n, ei and e.2 as shown in Fig. [I] The vector n is a constant unit vector pointing 
in the direction of angular momentum (1 = In). The two unit vectors ei and £2 in 
the orbital plane are oriented so that ei points in the direction of initial periapsis, 
apoapsis or toward the infinity (The choice depends on the orbit type and will be 




(1) 
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explained further in the text.). The components of these vectors with respect to 
the local Cartesian coordinate basis are expressed as in [21] : 

(cos u cos 1? — cos l sin w sin Q\ 
cos uj sin Q + cos i sin u> cos Q (2a) 
sin l sin to J 

(— sin uj cos Q — cos l cos o; sin ]?\ 
— sin oj sin Q + cos t cos uj cos (2b) 
sin l cos a; / 

h = (sin (.sin 12, — sin t cos ^2, cos t) (2c) 

where i7 is the longitude of the ascending node and t is the inclination of the orbit 
with respect to the X — Y plane (see Fig.^). 
By introducing a dimensionless variable 

u=— (3) 

and two dimensionless constants of motion related to orbital energy and orbital 
angular momentum |18j : 

2ME • u ^zjf 2M Y f ^ 

a = — - — and b = 2H I — — j (4) 

one can derive the differential orbit equation: 
du 



, . ±v'a 2 -w 2 (l- U ) + 6(l-u) , (5) 
dA 

where A is the true anomaly. As functions of u, time and proper time obey the 
following differential equations: 

dt 2Ma ,„ \ 

(6a) 



du u 2 (1 - u)^o? - u 2 (1 - u) + 6(1 - u) 
dr 2Ma 1 



du E v?y/a*-u 2 {\-u)+b(\-v) 



(6b) 
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Fig. 2 The polynomial P(u) and distribution of its roots in the interval u E [0, 1]. Orbits exist 
only where P(u) > 0, shown in colours. Corresponding orbit types are marked with letters A, 
B, C or D and the roots are marked with ui, 112 and M3. The sign of the discriminant D is 
also noted. 



After |5| is solved for it as a function of A, the orbit equation is written in vector 
form as 

2M 

r ( A) = — — (ei cos A + e 2 sin A) . (7) 

U(A) 

Solutions depend on the type of orbit, e.g. closed, scattering or plunging, and 
in the following section we present them for all types of time-like geodesicsrl 



2.1 Types of orbits 

Marking the polynomial in (JHJ with P(u) = a 2 — u 2 { \ — u) +6(1 — u), the solutions 
to |5]) — ( |6b[ ) exist only on intervals where P(u) > 0. This polynomial has three 
roots, while the discriminant D, which is defined as: 

a = 1 - 96 - y a 2 (8) 
/3 = -l-36 (9) 
D = a 2 +!3 3 , (10) 

determines the nature of these roots (i.e. the number of real/complex roots). Since 
orbits extend at most from u = to u = 1, only roots on this interval are of 
interest. In Fig. [2] the polynomial P(u) is plotted for all the four possible orbit 
types (according to the number of roots in the interval u G (0, 1)). 

The classification of orbits is more intuitive when it is done with respect to the 
effective potential V defined as |22j 

V =\j{l-u)(l + l 2 u 2 ) , (11) 

1 The differential equations (JEJl and | |6a[ | are formally the same for light-like geodesies |17| 
(for light-like geodesies take 6 = 0). 
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Fig. 3 The effective pot ential V for 
time-like geodesies (Eq. |ll| for 1 = 2.2. 
By choosing appropriate value of E, or- 
bits of any type can be constructed: or- 
bits with E > Vmax (blue area) are 
of type B (plunging orbits), orbits with 
1 < E < Vmax (yellow area) are either 
of type A (scattering orbits) or of type C 
(near orbits), orbits with V m in < E < 1 
(green area) are either of type D (bound 
orbits) or of type C, and orbits with 
E < V m i„ (red area) are only of type 
C. Note that for V m in < E < Vmax the 
discriminant is D < 0, and D > other- 



where I = 1/2M is the reduced angular momentum. Unlike in Keplerian case, the 
effective potential gains a maximum Vmax 



P \ (. 1 



' 1 + (P - Wl 2 - 3) 2 J V V ~ 



1/2 



at radius r ma Jf] 



2Ml[l-\jP -3 . (13) 



Clearly, the maximum exists only for I > \/3. The existence of Vmax greatly 
affects the nature of orbitsj^] especially if the orbital energy E is E ~ Vmax- such 
orbits can wind around the black hole at r max several times before continuing 
either away from or towards the black hole, and do not exist in case of Newtonian 
potential. If I = \/3, the maximum (and the minimum) of the potential disappears 
at r m ax = 6M, which is the radius of the last stable circular orbit. 

The effective potential and corresponding orbit types are shown in Fig. [3] The 
four types of orbits for massive particles have the following properties: 

- type A: scattering orbits with both endpoints at infinity. Scattering orbits can 
never extend below r = 3M. 

- type B: plunging orbits with one end at infinity and the other behind the 
horizon, 

- type C: near orbits with both ends behind the horizon of the black hole. 

- type D: bound orbits. Highly eccentric orbits can never reach below r = AM 
while circular orbits can never reach below r = 6M0 

Some typical examples of all types are shown in Fig. [4] Note that, only if E « Vmax, 
then r ma x corresponds to the radius of periapsis for type A and D orbits, and 
apoapsis for type C orbits. 



2 Note that r max is not the maximal radius an orbit can extend to, but the radius where 

Vmax — V {v-maxj- 

3 Obviously, the orbits can exist only for E > V . 

4 Highly eccentric orbits are orbi ts w ith energy E < 1 (which makes the orbits almost 
parabolic). From equations ( |12| and l |13| it follows that for type D orbits, r max is the smallest 
if Vmax ~ E ~ 1, which happens for I = 2 at r ma x ~ 4M. In this case, r ma x corresponds to 
the periapsis distance. 
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Fig. 4 Time-like geodesies with I = 2.32379 (for types A, B, and C) and I = 3 (for ype D). 
From left to right: Orbits of type A with E £ {1.0001, 1.035, 1.06, 1.083}, orbits of type B with 
E e {1.0887, 1.2, 1.6, 2.5}, orbits of type C with E £ {0.7, 0.97372899, 1.05, 1.086}, an orbit of 
type D with E = 0.988. The radius of the black circle is the Schwarzschild radius. 
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Fig. 5 The effective pot ential V for 
time-like geodesies (Eq.lTTJ for I = 1.9. 
By choosing appropriate value of E, only 
orbits of type B, C, and D can be con- 
structed: orbits with E > 1 (blue area) 
are of type B (plunging orbits), orbits 
with V m ax < E < 1 (top red area) 
are of type C (near orbits), orbits with 
Vmin < E < V m ax (green area) are ei- 
ther of type D (bound orbits) or of type 
C, and orbits with E < V m in (bottom 
red area) are of type C. Note that for 
Vmin < E < Vmax the discriminant is 
D < 0, and D > otherwise. 




Fig. 6 Orbit of type C for I = 1.9 and E = 0.99. The radius 
of the black circle is the Schwarzschild radius. 



If, however, I < 2 then V ma x < 1 (see (12)) and consequently, orbits of type 
A no longer exist, as shown in Fig. [5] Moreover, while orbits of type C still have 
both endpoints behind the horizon of the black hole, they can extend to infinity 
for E — > 1. An example of such extended type C orbit for I = 1.9 is in Fig.|6j 
Furthermore, if / is lowered below y/3, also type D orbits no longer exist, and only 
orbits of type B and C remain. 

Radial and circular orbits can be considered as special cases of type B and D 
orbits, respectively. The corresponding equations and parameters for radial orbits 
are: i = £7/(1 — u), r = y— (1 — u) + E 2 with zero angular momentum 1 = 0, while 



for circular orbits they are: r 
energy E = V mln (l) where V n . 



0, A = l/r = const., t = E/(l — u) = const. , wit h 



t (l) is the minimum of the effective potential (111 



2.2 Analytical solutions 

The solutions of the equations ^ - ( |6b[ ) are the following. 
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2.2.1 Types A and D 



In this case, the polynomial P(u) has either two (ui and 1*2) or three (ui, 712, 
and 113) real roots on the interval (0, 1) which can be elegantly expressed with 
the constants a, /3, and D using Cardan's formula [23] by introducing two more 
intermediary constants T> and ip [18] : 



ip = 2 arctan 



-D 



In terms of these, the roots can be written in the trigonometric form: 



til 



"2 



1 



l + 2|D|cos 
l + 2|D|cos 



u 3 = ■= 1 + 212)1 cos 



t 
3- 

t/} — 270 

3 / 

+ 270 



(14) 
(15) 

(16a) 
(16b) 
(16c) 



These roots can be associated to the radius of periapsis r p = 2M/u2 (types A, D) 
and cipOcipSlS Ta = 2M/U3 (type D only). Since the argument of arctan in ( |15[ ) is 
positive, it follows that < ip < tt, therefore u\ > U2 > 113. 
Using the substitution [23] 

u(x) = u 2 - (u 2 - U3) cos 2 x , (17) 

equations ([5| - ( |6b[ ) are transformed into Legendre form of elliptic integrals and 
integrated to obtain orbital variables A, t, and r as functions of x- 



A( X )=n(F(xH-K(m)) 



(18) 



t(x) 



2na 



I + U3 + 



n\ — m 



+ 



2(m — ni)(ni — 1) 

^ (E( X lm) - (l - -)F( X |m) 

(m — m)(ni — 1) \ V ni/ 

11 sin 2x\/ 1 — m sin 2 



77(m; x|m) + — ^— 77(n 2 ; xj™) 
1 — U3 



where: 



2(1 — n 1 sin 2 x) 
2n 



E * Z713 



iT(ni;x|m) + 



"3 



1 — 7*3 



n(n 2 ;x\m] 



m = 


7(2 — 713 
Ul - U3 


n = 


2 


\Ju\ — U3 


m = 


773 


712 = 


7i 2 - 7i 3 
1 — «3 



(19) 
(20) 

(21a) 
(21b) 
(21c) 
(21d) 
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r [M] 

Fig. 7 Left: Orbit of type A with E = 1.01 and I = 2.2. Black and red dots correspond to 
points at time intervals At = 5M and proper time intervals At = 5M, respectively. The black 
circle represents the Schwarzschild radius. Right: Time (black) and proper time (red) as a 
function of coordinate r, measured from the initial point at r = 34M. 



Inverting (181 by x(^) = am(K(m) + X/n \m) and using (171 one can also write 
the solution to the orbit equation ^ as a function of true anomaly in the form: 



u(\) = u 2 - (u 2 - u 3 )cn 2 (K(m) + -\m) . (22) 

n 

For type D orbits, both A and x can go from — oo to +oo. For type A, the values 
of X are in the interval \ £ (xmini 'Xmax')^ where \rain 

= arccos(^/u 2 /(u2 — uz)) 

and Xmax = arccos(— sju^jiui — 113)), while the values of A are in the interval 
A/n G (F ' (Xmin\iTi) — K(m), F(xmax\m) — K(m)). The values of x at periapsis and 
apoapsis are 7r/2 and respectively, while A = at periapsis. Definitions of elliptic 
integrals and functions are from Wolfram [25] . 

In Figures [7] and [HJ an example of solutions for a type A and type D orbits 
are shown, with black and red dots marking equal time and proper time intervals 
At = At = 5M. As expected, the dots are more widely spaced when closer to the 
black hole, and the lengths of sections corresponding to proper time intervals At 
are longer than those corresponding to time intervals At (which is also clear in 
the t = t{r) and r = r(r) plots of Fig. [7] and |8j where t{r) > r(r) for all r). Note 
that in case of type D orbit, the dots are plotted only for one orbital period, while 
the orbit is plotted for 3 periods to show the periapsis precession. 

In order to compare the efficiency and accuracy of the analytical expression 



( 19 1 for t to a direct numerical integration, equation ( 6a I has been integrated using 



fourth-order Runge-Kutta method with adaptive step-size control [26] . The elliptic 



integrals in equation (191 were calculated by Carlson's algorithm [27], while the 



Jacobi elliptic functions in (181 were from [26] . 

For type A orbit (E = 1.01, I = 2.2) the integration limits were r m 
6.15313M « r p and r max = 50M, while for type D orbits (E = 0.9704, 1 = I. 
the limits were r m i n = 5.04581M « r p and r max = 25.436M ~ r a . In both cases, 
the numerical integration fails, if r gets too close to either r p or r a since these are 



the zeroes of the polynomial P (u) in (6a|. Taking e.g. rmin = r p (l + 10 ) and 
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r [M] 

Fig. 8 Left: Orbit of type D with E = 0.9704 and f = 1.888. Black and red dots correspond 
to points at time intervals At = 5M and proper time intervals At = 5M, respectively. The 
black circle represents the Schwarzschild radius. Right: Time (black) and proper time (red) as 



a function of coordinate r, measured from the initial point at r ■ 



5.045M. 



rmax = f* a (l ~~ 10 -8 ) and thus avoiding the divergence]^] it turns out that numerical 
integration is ~ 50 — 80 times slower than analytical solution ( 19 ). In addition, the 
relative error 5t/t is ~ 4 orders of magnitude and ~ 2 orders of magnitude larger 



for numerical integration than for analytical solution ( 19 1 in case of type A and 
type D orbits, respectively. 

It should be also noted that some additional effort is required when numerically 



integrating (6a): if the orbit passes either r a or r p , e.g. a type D orbit spans many 



periods (or even just one!), or a type A orbit passes the periapsis, some book- 
keeping of periapsis and apoapsis passages has to be done in order to obtain 
the correct solution, e.g. by adding the correct number of half-periods. If using 



analytical solution, no such additional work is necessary, since ( 19 1 is essentially 
expressed with an angle along the orbit. 



2.2.2 Type B 



The polynomial P(u) can be factorized as P(u) = (u — ui)(u 2 + pu + q), where 
ui < is the only real root (see Fig.pl). The coefficients p, q, and the root ui are 



5 Consequently, the periapsis and apoapsis passage times have to be calculated in a different 
manner. 
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expressed as [23,18 



D : 

2) ■■ 

Ul ■■ 

p - 
q - 



(a-VD) 1 / 3 

b 

Ul — 1 

— b + pui 



: (l + D + D) 



(23a) 
(23b) 

(23c) 

(23d) 
(23e) 



Using the substitution [24 



«(x) = ui + \fu!{ 



2 X 

+ pui + q tan — , 



(24) 



equations ^ - ( |6b[ ) are transformed into Legendre form of elliptic integrals and 
integrated to obtain orbital variables A, t, and r as functions of X- 



A(x)=n(F(x|m)-F(x«,|m)) 



*(X) = 2a 



r(x) 



(25) 



+ 



+ 



(m -!)(! + fci) 
2\/\ni - m| 

2a 2 (m - 1) ( 1+ — 



(.♦ 




ni(7i 




fcl") 


Ql 




fcl " 









1 - ami 1 + 



1 



1 + 



1 — 7TI 



711—777 

F(x\m) 



Ql 



fcl 

In \xi 



n(m;x\m) 



77i sinxV 1 — 777 sin" x 
fci 1 — 771 sin 2 x 



E(x|m) - ^ ""^ v 1 A (1 + fci cosx) 



+ a 3 (l- 77 2 ) ( 1 + 



77(77 2 ;x|™) + 



A'2 



In 1x2 1 



(26) 



202 



2fci 



77l(l + fcl) / fcl 771 + 1 - 2777 

H ,. I ^ 7 2 In si + 



2 \/l ni ~ m l ~ 1 m-m 

2 



fcl - 1 



F(x|m) 



+ 



1 



777 — 771 V fcl — 1 



fcl 



E(xl' 



77i smx 



fci 1 — 771 sin x 



—2— (1 + fci cosx) 



(27) 
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where: 



(Ui +pui + q) 
m 3 



ui + p/2 
\fu\+ pu\ + q 
1/4 



Oil 



XI 



X2 



n 2 ui + 1 



«2 



1 



1 + u\n? 
k\ - 1 



ri2 



(n 2 ui + l) 2 
1 + (1 - Ul )n 2 
1 - (1 - ui)n 2 

kl 



a 3 



n 2 (l — ui) — 1 



k ; 2 



m sin x 



m sin x 



y/n2~ 



m Bin x 
m sin x + \Zl~~ 



m sin x 



m sin x 



(28a) 

(28b) 
(28c) 

(28d) 
(28e) 
(28f) 
(28g) 



m sin x 

Inverting (251 by xW = am (F(X' 



m sm x 



\m) + X/n \m) and using (24) it is straightfor 
ward to obtain the following form of the orbit equation: 

1 1 - cn(F(xoo|m) + ±\m) 



t(A) = ui + 



l + cn(F( X oo|m) + ^|m) 



(29) 



The values of x are m the interval x € {xbh, X°o)> where X-B/f = arccos 



l-n 2 (l- 
l+ra 2 (l- 



and Xc 



-n 2 u 1 



Since neither periapsis nor apoapsis exist for this type 



of orbits, the value of A is measured from the direction toward infinity, i.e. A = 
atr^oo and the values of A are in the interval X/n £ {F(xBH\m) — F(xoo|m),0). 



Additionally, it is clear from equations (6a) and (6b| that while time t diverges as 



r —> 2M , proper time r remains finite (see Fig.[9). 

In Fig.|9j an example of the solution for a type B orbit is shown, with black 
and red dots marking equal time and proper time intervals At = At = 1M . As in 
previous case, the dots are more widely spaced when closer to the black hole, and 
the lengths of sections corresponding to proper time intervals At are longer than 
those corresponding to time intervals At. However, since t — > oo when r — > 2M, 
the black dots start to concentrate at r ~ 2M, while the red ones remain distinctly 
separated. This is also visible in the t = t(r) and r = r(r) plots of Fig.[9j where 
t(~ 2M) diverges and r(~ 2M) has a finite value. 



The efficiency and accuracy of the analytical expression ( 26 1 for t compared 



to a direct numerical integration of (6a) has been done using the same methods 
as in the previous case. For type B orbit (E = 1.06, I = 2.2) the integration limits 
were r m i n = 2.0001M and r max = WOM. The numerical integration is ~ 20 times 
slower than analytical solution (26) and the relative error 5t/t is ~ 2 orders of 



magnitude larger for numerical integration than for analytical solution ( 26 1 



2.2.3 Type C 

Since type C orbits exist for both D > and D < 0, two different sets of pa- 



rameters are introduced: If D > 0, use the parameters (23) and (28) for type B 
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Fig. 9 Left: Orbit of type B with E = 1.06 and I = 2.2. Black and red dots correspond to 
points at time intervals At = 2M and proper time intervals At = 2M, respectively. The black 
circle represents the Schwarzschild radius. Right: Time (black) and proper time (red) as a 
function of coordinate r, measured from the initial point at r = 29M. 



orbits. If D < 0, use the parameters from ( 14) - ( |16[) to calculate p = —(112 + 1*3) 
and q = U2U3, and use them in (23d I - (23e| and (|28[). In both cases, the root ui 
can be associated to the radius of = 2M/u±. Also, if m > m, do the 



following substitution in equations ( 26 1 and ( 27 ) 



In |:ei I — > 2arctan(yi) 



where y\ is 



yi 



n\ sin x 



(30) 



(31) 



m sin x 

This substitution is necessary because if m > n\ then x\ becomes complex, so 
it is more convenient to use the relation ln((l + ix)/(l — ix)) = 2Artanh(ia;) = 
2i arctan(a;), where the imaginary unit i cancels out with i from y/m — n\ in front 
of the In term. 

While the solutions for u, t, and r are the same as for type B, the solution for 



A is 



Kx) = n¥{x\m) 



(32) 



i.e. use equations ( 24 ) - ( 29 ) with the above replacements. Inverting ( 32 1 by xW = 
am(A/n |m) and using (24 1 it is straightforward to obtain the following form of 
the orbit equation: 

' A I 



i(A) = u\ + 



1 1 



:\ m ) 



l + cn(£|m) 



(33) 



l+n A (1 - 

In case of 



The limits for A are X/n 6 {— F(xBH\m), F(xbh|to)}, where xbh = 
Note that in this case, the values of A and x a t apoapsis are A = x = 
type C orbits it is also true that for r — > 2M, time t diverges and proper time r 
remains finite. 
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Fig. 10 Left: Orbit of type B with E = 1.1 and I = 2.8. Black and red dots correspond to 
points at time intervals At = 0AM and proper time intervals At = 0AM, respectively. The 
black circle represents the Schwarzschild radius. Right: Time (black) and proper time (red) as 
a function of coordinate r, measured from the initial point at r = 2.0001M. 



In Fig.|10| an example of the solution for a type C orbit is shown, with black 
and red dots marking equal time and proper time intervals At = At = AM. 
As in case of type B orbit, since t — > oo when r — > 2M, the black dots start to 
concentrate at r ~ 2M, while the red ones remain distinctly separated. This is also 
visible in the t = t(r) and r = r(r) plots of Fig. [9] where t(~ 2M) diverges and 
t(~ 2M) has a finite value. Note that since the orbit is always very close to the 
black hole, the difference between t and r is huge, so the number of At intervals 
is much smaller than the number of At intervals. 

The efficiency and accuracy of the analytical expression ( 26 1 for type C orbits 
compared to a direct numerical integration of (6a) has been done using the same 
methods as in previous cases. For type C orbit (E = 1.1, I = 2.8) the integration 
limits were r m i n = 2.0001M and r max = 2.50581839M « r a . As in the case of type 
A and D orbits, the numerical integration fails, if r gets too close to r a . Taking 
e.g. rmax = r a (l — 10 -8 ) to avoid the divergence, it turns out that numerical 
integration is ~ 270 times slower than analytical solution (261 and the relative 
error St/t is ~ 2 orders of magnitude larger for numerical integration than for 
analytical solution (26 I for type C orbit. Also, similarly as in case of type A and D 
orbits, if the orbit passes the apoapsis r a , this has to be done taken into account 



only if doing numerical integration of ( 6a 



3 Summary 

In this paper, the analytical solutions of the orbit equation for time-like geodesies in 
Schwarzschild space-time are presented in a very straightforward way. The orbits 
are classified into four types according to the roots of polynomial P(u). This 
classification is also presented in a more intuitive way, i.e. according to the effective 
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potential and orbital energy. The four orbit types are: type A - scattering orbits 
with both endpoints at infinity, type B - plunging orbits with one end at infinity 
and the other behind the horizon, type C - near orbits with both ends behind 
the horizon of the black hole, and type D - bound orbits. The analytical solutions 
are expressed with Jacobi elliptic functions where the true anomaly is the only 
parameter. 

The analytical solutions for time and proper time for all four orbit types are 
also presented here and are expressed as functions of one parameter %. A simple 
relation between x an d true anomaly A is given for all four types. 

Since these analytical solutions for time and proper time are expressed with 
elliptic integrals, which can be numerically calculated very efficiently and accu- 
rately either with Landen transformations [23] or Carlson's algorithms [27], they 
can be very useful in particular for modelling dynamical phenomena near black 
holes. These solutions have been in fact already successfully used together with 
light-like solutions 17 in modelling tidal disruption of low-mass satellites around 
black holes [28] and quasi-periodic oscillations from X-ray binaries [29] . Although 
the motivation for this work comes from black hole physics, the method was se- 
lected due to its performance [13] also for investigation of a relativistic approach 
to Galileo Global Navigation Satellite System [21] . 
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